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^ Abstract 

We adopt an operatorial method based on the so-called creation, annihilation and 
— number operators in the description of different systems in which two populations 

• • 

^ interact and move in a two-dimensional region. In particular, we discuss diffu- 

sion processes modeled by a quadratic hamiltonian. This general procedure will 
be adopted, in particular, in the description of migration phenomena. With re- 
spect to our previous analogous results, we use here fermionic operators since they 
automatically implement an upper bound for the population densities. 
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1 Introduction and preliminaries 



A large body of theoretical and experimental evidence that spatial patchy environments 
influence the dynamics of species interactions is available in the literature [TJ [2j [3j HJ El E] . 
Hence, a broad variety of spatially extended models has been developed in theoretical 
biology. The principle of competitive exclusion (Gause's Law), stating that two like 
(identical) species cannot coexist in the same habitat, is violated in patchy environments 
where two like species may coexist due to migration [7]. A lot of evidence exists about 
the coexistence, as met apopulat ions jH] in a fragmented environment, of two competing 
species (or populations) even if one is competitively superior to the other; in fact, the local 
extinction in a patch of the inferior competitor has no global effect if it is able to disperse 
more effectively into unoccupied patches (see [9] and the references therein quoted). 

Besides the usual models based on continuous reaction-diffusion equations [TU] and cel- 
lular automata, the coupled map lattice (CML) formalism has been widely used in the sim- 
ulation of biological spatial interactions. In the usual CML approach, local prey-predator 
(or host-parassitoid) dynamics are coupled with their n-nearest neighbors through some 
appropriate exchange rule. Populations interact and disperse over the points of a lattice 
(used to simulate the patchy environment). In such a context various aspects can be 
observed, such as the emergence of some persistent spatial patterns in the distributions of 
the competing species (e.g., phytoplankton distribution in the oceans), or the phenomena 
of synchronization between the phases of nearby regions [TTj . 

The mathematical oriented literature on migration is very poor. But for those cited 
above (more concerned with the coexistence of biological species), only few papers deal 
with this problem. For instance, in |12j . the authors propose a game-theoretic model of 
migration on animals, while in [13] . the author concentrates his attention to migration 
in The Netherlands. Moreover, it is worth of note the paper [13] where an oscillating 
behavior has been described with Sheppard's generalized dynamic migration model. 

Other models of competing populations including also spatial interactions can be con- 
structed using completely different tools, say, operator methods of quantum mechanics. 
In a series of recent papers one of us (F.B.) used some tools from quantum mechanics, like 
operator algebras and, in particular, the so-called number representation, to discuss some 
toy models of stock markets pj3 EH [T71 [18] . More recently, we have adopted the same 
framework in the analysis of love affairs, [HI [201 EI] ■ In these rather different areas the use 
of the number representation proved to be quite useful to set up a natural description of 
the system under consideration. In fact, several quantities which change discontinuously 
are very well described in terms of the integer eigenvalues of certain relevant self-adjoint 
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operators, the observables of the system, and the dynamics is described by an energy-like 
operator, the hamiltonian. This approach has produced several interesting results like, 
just to cite a few, the possibility of describing the time evolution of the portfolio of the 
traders of a simplified stock markets, or of computing the transition probability of the 
market itself from a given initial state (i.e., from a given distribution of the cash and 
the shares for the different traders) to a final state. In the context of love affairs, we 
have found a non-trivial dynamical behavior of the actors of the affairs (Alice and Bob), 
even in the presence of a third actress, Carla, and we have also considered the effect of a 
reservoir interacting with them, mimicking the real world. 

In this paper, we use an analogous strategy, based on fermionic rather than bosonic 
operators, to describe some models arising in different contexts. In particular, we will 
describe a (strictly local) situation in which two populations live together and are forced 
to interact, and a second situation in which the two species occupy (in general) different 
cells of a two-dimensional lattice, interact and move along the cells. These models can 
be useful to model many different systems. In particular, we will restrict ourselves to 
consider the following two: a migration process in which a population moves from a given 
(poor) place to a richer region of the lattice, which is already occupied by a second group 
of people, and a system in which the two populations compete. 

The main motivations which suggest the use of the tools originally developed in a 
quantum context to describe classical situations have been widely discussed along the 
years, [151 EH EHl HI HH I2Q1 121], and will not be repeated here. We just want to mention 
that, in recent years, the intersection between quantum and classical worlds has became 
more and more frequent. We refer to [231 1211 125] , and to the references therein, for many 
other applications and for some general comments. 

The choice of using fermionic operators rather than bosonic ones is mainly based on 
two reasons. The first one is of technical nature: the Hilbert space of our models is 
automatically finite dimensional, so that all the observables are bounded operators. The 
second reason is related to the biological and/or sociological interpretation of our model: 
for each population which will be considered, we will have two only possible non-trivial 
situations. In the first one (the ground state) there is a very low density, while in the 
second one (the excited state), the density is very high. Hence, if we try to increase the 
density of the excited state, or if we try to decrease the density of the ground state, we 
simply annihilate that population! We can interpret this fact just saying that there exists 
upper and lower bounds to the densities of the populations which can not be overcame 
for obvious reasons: for instance, because the environment can not give enough food 
to the populations. Of course, this rather sharp division in just two levels may appear 
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unsatisfactory. However, it is not hard to extend our procedure to an arbitrary number of 
levels, paying the price of some technical difficulties. We will not consider this extension 
here since, already in our hypotheses, an absolutely non-trivial and realistic dynamics 
can be still deduced. 

The paper is organized as follows. In Section [2j we consider a first simple model 
involving two populations, and we analyze the dynamics of their relationship starting 
from very natural assumptions. The considered model is linear, strictly local, and the 
equations of motion are solved analytically. 

In Section |3j we extend this model by allowing a spatial distribution. The interaction 
is quadratic, so that the solution can be again deduced analytically. The model will be 
discussed in terms of migrant and resident populations. 

In Section [4], we show how the same model introduced in Section |3j with a different 
choice of the parameters and of the initial conditions, can be used in the description of 
two competing populations. 

Finally, Section [5] contains our conclusions, while, to keep the paper self-contained, 
we list in the Appendix few useful facts about quantum mechanics and the number repre- 
sentation for fermions. In all the considered cases a quasi-periodic dynamics is obtained 
and the two populations coexist in the patchy environment. 

2 A first model 

In this section, we introduce a first simple model which is useful to fix the main ideas and 
the notation. This model will also be used in the next section as a building block of a more 
sophisticated model. In particular, no spatial distribution will be considered. Following 
[T5l IT6j \T7\ \18\ [191 [20l |2T] , we associate to each population Sj of our model an annihilation 
and a creation operator aj and aj, and a related number operator hj := aja^. Here, we 
just consider two populations, Si and i?2- We assume the following anticommutation 
rules: 

{a h a]} = 5 id , {a h a,} = {at, a]} = 0, (2.1) 
i,j = 1,2. Recall that {x,y} = xy + yx. These rules imply in particular that a 2 - = 

(a^j = (see the Appendix). Hence, if <^o,o is the ground state, aiv?o,o = 02^0,0 = 0, the 
only non-trivial vectors of our Hilbert space % are 

^0,0, y?i,o := 4^o,o, <A),i := 4<A),o, ¥>i,i : = «i«2^o,o- 
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This means that dim(H) = 4. The biological interpretation of these vectors follows from 
the following eigenvalue equations: 

ftlVnum = fWni.rw,) ™2<Prn,n 2 = ^^rii.nj, ( 2 -2) 



j = 1,2, where = ctjCij is the number operator of Sj (see the Appendix). Saying 
that </9o,o is ^e s^^e 0/ t/ie system means that there are very few subjects of the two 
populations in our region. If the state is <^i,o, then there are very few elements of S 2 but 
very many elements of «Si. The opposite situation is described by </? 0) i, while ipx,i describes 
the case in which both populations are abundant. As already stated, it is not possible 
to have, e.g., more elements of Si than those described by ipi$ or (fix- trying to further 
increase the density of Si simply destroys this population! This is a simple consequence 
of (aj) 2 = 0. As in our previous applications, we use a self-adjoint operator, which we call 
the hamiltonian of the system, to derive the dynamics of the system, and which describes 
the interaction between the populations. The self-adjoint hamiltonian which we assume 
here is the following one: 

H = H Q + \H I , H = Uia\ai + co 2 a 2 a 2 , Hi = a\a 2 + a 2 a i> (2.3) 

in which uj and A are real positive quantities. In particular, it is A = when the two 
populations do not interact. In this case, H describes a static situation, in which the 
densities of the two populations, described by the number operators hj, do not change 
with t. This is a consequence of the fact that [H , hj\ = 0, j = 1, 2. On the other hand, 
if A 7^ 0, then Hi describes a situation in which the density of Si increases (because of 
a\) while that of S2 (decreases because of a 2 ), or viceversa, [22]. The equations of motion 
that are obtained are (see (A. 2) in the Appendix): 

ai(t) = —iuiai(t) — i\a 2 (t), ^ 

a 2 {t) = -2w 2 a 2 (t) -zAai(t), 

which can be solved with the initial conditions ai(0) = ai and a 2 (0) = a 2 . The solution 
looks like 

ai(t) = 777 (ai ((wi - w 2 )$_(*) + ^+(t)) + 2Aa 2 $_(*)) , 

2 / (2.5) 

a 2 {t) = y 5 ( fl 2 (-K " w 2 )M*) + + 2Aai$_(*)) , 

where 



5 = v /(u; 1 -w 2 ) 2 + 4A 2 , 

= 2exp f 1 ^ 2; j cos ^- 

$_(*) = -2zexp ( ^ ^ J sin ( - 
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It is now easy to deduce the mean value of the time evolution of the number opera- 
tor hj(t), which, as discussed before, we interpret here as the density of Sf nj(t) := 
({p nit n 2 ,hj(t){p nitn2 ). More explicitly, rij{t) is the time evolution of the density of Sj as- 
suming that, at t — 0, the density of Si was n\ and that of S 2 was n 2 , the quantum 
numbers labeling the state f ni ,n 2 - Using (2.2) and the orthonormality of the different 
V^ni.^'s, we obtain 



7ll(t) = Til 

and 
n 2 {t) = n 2 



lux - u 2 ) 



+ 



4A 2 



{u x - uj 2 ) 2 + 4A 2 (ut - uj 2 ) 2 + 4A 2 
(wi - w 2 ) 2 



rii cos 



+ 



4A 2 



(wi - w 2 ) 2 + 4A 2 (wj - w 2 ) 2 + 4A 



^2 cos 



n 2 sin 



ni sin 



(2.6) 



(2.7) 



Notice that these formulas imply that ni(t) + n 2 (t) = n\ + n 2 , independently of t and A. 
This is expected, since it is easy to check that [H, hi + h 2 ] = 0. Secondly, since n\ and n 2 
can only be or 1, we can also check that, if n\ = n 2 = n, then n\(i) = n 2 (t) = n for all 
t: if the two populations are equally distributed at t — 0, then they do not change their 
distributions. If % = 1 and n 2 = 0, then 



, 4A 2 . 2 
1 =- sin 



while if n\ = and 712 = 1, then 

n 2 (t) 



, 4A 2 . 2 
1 ^r- sin 

0^ 



n 2 (t) 



m(t) 



4A 2 
<5 2 



4A 2 
5 2 



sin 



sin 



In all these cases we have < rij(t) < 1 for all t, as it should be. 

Since n\(t) + n 2 (t) = n\ + n 2 , we find that Sn(t) := \ni(t) — rii\ 
in the two cases above, gives 

4A 2 



\n 2 (t) — n 2 \ which, 



Sn(t) 



sm 



so that the variations of the two populations coincide. In general, equation (2.6) gives 

4A 2 



ni(t) - m 



{n 2 — rii) sin 



which is in agreement with the previous result since n\ and n 2 can only be or 1. Notice 
also that, if ri\ = n 2 , then ni(t) = n\ for all t, as already remarked. 
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Now let us restrict ourselves, for concreteness, to the case ri\ = 1 and n 2 = 0. Hence 

r c / M 1 r t(5 7T 37T 

A n := maxj onit) \ = ^, tor — = — , — , 

1 WJ 1 _|_ { Ul ~^ 2 ) 2 2 2 

while 

min{5n(t)} = for — = 0, 7r, 27r, .... 

In particular, A n is almost equal to 1 if oj\ ~ u 2 , independently of A ^ 0, while is 
almost 0, as A is kept fixed, when \u\ — u 2 \ is very large. Incidentally, if A = the two 
species do not interact and, in fact, A n = 0: the model has essentially no dynamics. For 
obvious reasons we are only interested to the case A > 0. As already mentioned before, 
these results show that the free hamiltonian, which does not affect the density of the 
populations if A = 0, produces a non-trivial effect if A 7^ 0. A similar conclusion will be 
deduced in Sections [3] and [4], where we consider the spatial version of this model. More in 
details, \oj% — uj 2 \ can be considered as a sort of inertia of the system: the larger its value 
is, the smaller the variations of rij{t) — rij are. On the other hand, if \u)\ — u 2 \ — 0, then 
the system has a very low inertia and, in fact, very large changes in the densities of both 
populations are allowed. An interesting fact is that only the difference between the two 
frequencies oj\ and uj 2 play a role in the dynamics of both Si and S 2 . 

Concerning a relation between the w's and A, from (2.6) and (2.7) we also see that if 
\oji — uj 2 \ ^> 2A there is essentially no dynamics: n,j(t) ~ rij, j = 1,2. On the contrary, 
if \uj\ — W2I 2A, the constant contributions in (2.6) and (2.7) are very small compared 
with the oscillating contributions. These results agree with our previous conclusions. 



3 A spatial model 

In this section, we extend the model introduced above with the aim of including spatial 
effects: we consider a 2D-region 1Z in which, in principle, the two populations are dis- 
tributed. Under reasonable assumptions, a simple model for <Si and S 2 can be deduced, 
and its dynamics investigated. 

The starting point is the {e.g., rectangular or square) region 1Z, which we divide in N 
cells, labeled bya = 1,2, . . . , N = L ■ L' . With a = 1 we label the first cell down to the 
left, while iV is the last cell, up to the right (see Figure [T]). 

The main idea of our model here is that in each cell a the two populations, whose 
related operators are a a , o) a and = a^ a a a for what concerns <Si, and b a , tf a and ni 6 ' = 
b^ba for <$2, behave as in the previous section. This means that the same hamiltonian as 
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Figure 1: The two-dimensional lattice for the spatial model. 



in (2.3) will be assumed here in a. Using our new notation, we define 

H a = H° a + \ a Hi, H° a = u a a o) a a a + u b a blb a , Hi = a\b a + b{a a . 



(3.1) 



Extending what we have discussed in the previous section, it is natural to interpret the 
mean values of the operators and as local density operators (the local densities are 
in the sense of mixtures; hence, we may sum up local densities relative to different cells) of 
the two populations in the cell a: if the mean value of, say, n~a \ in the state of the system 
is equal to one, this means that the density of S\ in the cell a is very high. Notice that 
H a = iff, since all the parameters, which in general are assumed to be cell-depending 
(to allow for the description of an anisotropic situation), are real and positive numbers. 
The anticommutation rules are those in ( 2.1| , which we rewrite as 



{a«,aU = {b a ,bl} = 5, 



(3.2) 



Of course, the hamiltonian H must consist of a sum of all the different H a plus another 
contribution, h, responsible for the diffusion of the populations all around the lattice. A 
natural choice for h is the following one: 



h = ^2p a ,f3 1 la (a a a^ + apal) + 7 b (b a b^ + b p b ] a 



(3.3) 



where also j a , 76 and the p a ^ are real quantities. In particular, p a ^ can only be or 
1 depending on the possibility of the populations to move from cell a to cell (3 or vice- 



versa. In fact, this is the meaning of, e.g., the contribution a a at in (3.3), [22j: because 
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of a a the mean value of n£ decreases, while the mean value of increases because 
,t 



of dp. We interpret this feature saying that some individuals of Si are moving from the 
cell a to the cell (3, if p a $ = 1. For this reason the p Q ^s can be considered as diffusion 
coefficients. Notice that a similar role is also played by 7 a and 7&, which, however, could 
assume different values. This is important since, in this way, a different mobility for the 
two species can be introduced. In the rest of the paper we will assume that diffusion 
may take place only between nearest neighboring cells. Of course, we should clarify what 
we mean by neighboring: in this paper, we will consider a simple planar topology, in 
which the neighboring cells of the cell labeled a are the cells a — 1, a + 1, a + L and 
a — L, provided that they exist (this is verified for the internal cells of the lattice); for 
the cells located along the boundaries of the lattice we have only three neighbors (two 
neighbors for the four cells located at the corners of the lattice). This is a natural choice 
for the physical system we have in mind here. However, different choices could also be 
considered. For instance, we could also use a torus topology in which all the cells have 
four neighboring cells. To deal with these different topologies, it is enough to modify 
the diffusion coefficients. We will assume that p ajCe = and that, for symmetry reasons, 
Pa,p = Pp,a- Now we define H = J2 a H a + h, which is self-adjoint. The differential 



equations for the annihilation operators (see (A. 2) in the Appendix) read 

a a = -%u) a a a a - i\ a b a + 2z'7 ^Jp^g^a, 

* (3-4) 
b a = -iu b a b a - i\ a a a + 2ij b } j p a> pb /3 . 

a 

Remark 1 Readers with a background in quantum many-body could interpret the operator 
H as the hamiltonian of two kinds of fermions mutually interacting. From this point of 
view, it could be interesting consider other aspects of the system, and not only its dynamics. 
Phase transitions, Green's functions, non-zero temperature states are all typical problems 
which could be of a certain interest starting from a similar hamiltonian, J2i 



Remark 2 It might be worth noticing that the fermionic operators could be replaced by 
Pauli matrices. Using these operators (which again can be considered as raising, lowering 
and number operators) we would obtain a formally different hamiltonian describing the 
same physics. 

3.1 A simple case: equal coefficients 

As a first step, we suppose here that ui h a = = u, X a = X and r y a = lb = 7, for all a £lZ. 
Hence, by introducing a a (t) = A a (t)e~ zwt and b a {t) = B a (t)e~ tuJt , the above equations can 
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be easily rewritten as 

A a = -i\B a + 2ij22p atl3 Ap, 

' (3-5) 
B a = -i\B a + 2ry } ( p a ,pBp, 

P 

which are true independently of the size of the region 1Z. 

The diffusion coefficients p a ^ will all be zero but when a and refer to nearest 
neighbors in the planar topology. In this case p a> p = 1. We will consider now the 
situation of a square region 1Z with N = L 2 cells, starting with the simplest non-trivial 
situation, L = 2. In this case, the only non zero diffusion coefficients are pi >2 , Pi,3, P2,i, 
P2,4, P3,i, P3,4, P4,2 and p4 ; 3, all equal to one, while the remaining ones are zero. Fixing 
natural initial conditions, A a (0) = a a and B a (0) = b a , we get, for instance, 

Ai(t) = -[2(ax — a 4 ) cos(t) + (aj + a 4 + b 2 + b 3 ) cos(7t) + (a x + a 4 — b 2 — b 3 ) cos(9t) 
- i(2(6i - & 4 + (&i + 6 4 ) cos(8t)) sin(t) - (a 2 + a 3 )(sin(7t) + sin(9t)))], 

and so on. The number operators (i.e., the local densities of Si and S 2 ) are deduced di- 
rectly from the capital operators Aj(t) and B 1 (t), since n^\t) = ed(£)a 7 (i) = A^,(t)Ay(t) 
and hj (t) = 6^,(t)6 7 (t) = BUt)Bj(t). Assuming that for t = both the populations 
are concentrated in the cell I 1 , (0) = n[ (0) = 1, while n^(0) = ni^(0) = for 
a = 2,3,4, we find the following results 

n[ a \t) = nf\t) = (cos(4t)) 4 , 

n 2 a) (t) = 4 6) W = ^(sin(8t)) 2 , 

f (3.6) 
4 a) W = 4" ) W = i (sin(8t)) 2 , 

4 a) W = ^W = (sin(4t)) 4 . 
These results look reasonable because of the following considerations: 

1. taking into account the fact that, because of our simplifying assumptions, all the 
parameters of a and b coincide, and since their initial conditions coincide as well, it 
is clear that the spreading of the two populations must be identical; 

2. at t — only the cell 1 is populated; 



1 Of course this choice is not very relevant in the context of migration but is useful just to fix the ideas. 
The application to migration will be considered for a larger lattice, where the situation is surely more 
realistic. 
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3. we observe that n£\t) = n^\t) and nf\t) = nf\t); this is not surprising since, 
because of the isotropy of 71, there is an equal probability for, say, a member of Si 
to move from cell 1 to cell 2 or to cell 3; so he can reach cell 4 only through cells 2 
or 3, but not directly; 

4. this explains why, if we plot n[ (t), (t) and nf 1 if), we see that, for small values 
of t, n 2 if) increases faster than n± if) but, after some time, the opposite happens; 
this is because after cells 2 and 3 are populated, they both start contributing to the 
population of cell 4. 

Let us now move to L = 3. In this case TZ is made up of 9 cells, and the only non-zero 
diffusion coefficients are (listing just one between p a ^ and pp ta ) Px,2, Pi,4, Pz,5, P2,3, P3,e, 
Pi,5, P4,7 ; P5,6> P5,8) P6,9; P7,8 an d P8,9> which are all equal to 1. The 18 differential equations 



extending those in (3.5) can be written as 



Xq = iMqXq 



M 9 = 27M 9 - AJ 9 , 



(3.7) 



where we have introduced the following vector and matrices: 



Xq 



A 2 



Aq 

B 2 



\Bq ) 



Nq 



/oioiooooo\ 

1 1 1 

1 1 

1 1 1 
10 10 10 10 
1 1 1 
1 1 
1 1 1 
1 1 



Mq 



Nq Oq \ j = ( 9 1 9 

Oq N 9 )' \ h 9 

where 9 and 1 9 are respectively the 9x9 null and identity matrices. Notice that Aiq is 
a symmetric real matrix. 

The generalization to larger TZ is straightforward. In this case we have 



X L 2=iM L 2X L 2, M L 2 = 2^M L 2 - XJ L 2. 



(3.8) 
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M L 2 



Here the transpose of X L 2 is (Ax, A 2 , ■ ■ ■ , A L 2, B x , B 2 , . . . , B&), while L 2, 1 L 2 and J L 2 
extend those above, and 

' N L 2 L 2 
L 2 N L 2 

Once again, M.l 2 is a symmetric real matrix. Of course, the explicit form of the matrix 
Ntf can be constructed extending the previous considerations: this matrix have all zero 
entries but those matrix elements corresponding to nearest neighbors, which assume as 
values 1. 



The solution of equation (3.8) is 



X L 2 (t) = exp (i M L 2t) X L 2 (0) . 

Let us call d a ^(t) the generic entry of the matrix exp(z Ai i?t) , and let us assume that 
at t — the system is described by the vector y? n a n f,, where n a = (nj, n%, . . . ,n a L2 ) and 
n b = (n\, n\, . . . , n b L2 ) ■ Hence, the mean values of the time evolution of the number 
operators in the cell a, 

N a(t) = (^n»n» ) fl t a Wfl tt (% n »n*) = (W.nf , A a (t)A a (t)(p aa>I1 b ) , 
N a(t) = (<Pn'>,n'>,b f a ( t ) b cx(t)<Pn",n>>) = (<Pn",n>>, B l(t)B a (t)<p n a^ b ) , 

can be written as 



N a a (t) = \d a ,e(t)\ 2 n a e + £ \d a , L , +e (t)\ 

0=1 8=1 
L 2 L 2 

N b a (t) = Yl \d^+afi(t)\ 2 n a e + \dL 2 +a ,L 2 +e(t)\ 



(3.9) 



2 nl 



3.2 Back to the general case: migration 



The same strategy which produces solution (3.9) can be used to solve system (3.4). In 



this case, (3.8) is replaced by a similar equation, 

X L 2 = %K, L 2X L 2, (3.10) 



where /C^2 = 2T L 2 — P L 2, with T L 2 and Pip. two L 2 x L 2 matrices defined as follows: 

T L 2 = 



N L 2 with 1 replaced by 7 a 

N L 2 with 1 replaced by % 
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Figure 2: Evolution of local densities (solid line for Si and dashed line for 52). Africa: top-left (a); 
Europe: top-right (b); Mediterranean Sea: bottom (c). j a — 0.1, 7& = 0.004, w° — 1, uj h a = 0.3, 
\ a = 0.05, Va G K. 



local density 



local densily 





Figure 3: Evolution of local densities (solid line for S\ and dashed line for £2)- Africa: top-left (a); 
Europe: top-right (b); Mediterranean Sea: bottom (c). j a — 0.1, 76 = 0.004, w° = 1, uj b a = 1, X a = 0.05, 
Va e H. 
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Figure 4: Evolution of local densities (solid line for Si and dashed line for 6>2). Africa: top-left (a); 
Europe: top-right (b); Mediterranean Sea: bottom (c). j a — 0.1, 7b = 0.004, w° = 1, = 3, X a = 0.05, 
Va e K. 



and 



Pl 2 



fiW A 



a 

Here we have introduced the following matrices: = diag{w°, wf, • • • , w 2 2 )' ^ 

di&g{u\,uj%, . . .,u b L2 }, and A = diag{Ai, A 2 , . . . , \&}. 



The solution of equation (3.10) is 



X L 2 (t) = exp (i K L 2t) X L 2 (0) . 

Calling f a , pit) the generic entry of the matrix exp(i/C L 2t), and repeating the same pro- 
cedure as above, we get 

L 2 L 2 

K(t) = \Ue(t)\ 2 n a e + \UL 2 + e(t)\ 2 n b e , 

(3.11) 

N b a (t) = J2 \fLH«At)\ 2 n<e + E \fL> +a ,LZ + e(t)\ 2 n b e . 

0=1 9=1 

These formulas are used to deduce the local densities of the two populations S\ and S2 
in three different regions. The first one, 1Z\, corresponding to cells 1, 2, L + 1 and L + 2 
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Figure 5: Evolution of local densities (solid line for Si and dashed line for 52). Africa: top-left (a); 
Europe: top-right (b); Mediterranean Sea: bottom (c). j a = 0.1, 76 = 0.004, uj^ = 1, = 0.3, Va G 7Z; 
\ a = 0.2 for oe^URz and A Q = 0.05 for a € 7£ 3 . 

(bottom-left corner of 72.), is that part of TZ where all the members of Si are originally 
(i.e., at t = 0) localized. Population 1S2, at £ = 0, is assumed to be localized in the four 
cells L 2 — L — 1, L 2 — L, L 2 — 1 and L 2 , the region 72-2 (top-right corner of TZ). All the 
other cells form 72-3, that part of TZ which must be crossed by the populations to reach 
the other region of the lattice. Just to fix the ideas, we could think of Si and 52 as people 
from Africa (TZi) and Europe (72-2), respectively, and the Mediterranean sea as the region 
72.3. We also fix L = 11. In Figures [2]-[7] we plot the two local densities (the sum of the 
densities in the different cells) for «Si (solid line) and S 2 (dashed line), in TZi (top-left 
plot), TZ 2 (top-right plot) and TZ 3 (bottom plot), for different choices of the parameters 
and for the same initial conditions given above. 

Remark 3 The reason why we are talking here of Africa and Europe is very much related 
to what we have experienced in Italy, and in Sicily in particular, during this last year, with 
all the people moving from Africa, and from Libya in particular, and reaching Lampedusa 
first, and Europe soon after. In the past century, to a similar migration process took part 
a lot of people coming from Sicily ( among whom many of our relatives ) and moving to 
America looking for a better life. Many of them reached a reasonable well-being, and some 
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Figure 6: Evolution of local densities (solid line for Si and dashed line for 6>2)- Africa: top-left (a); 
Europe: top-right (b); Mediterranean Sea: bottom (c). -f a — 0.1, 7b = 0.004, w° = 1, w b a — 1 Va e 1Z; 
\ a = 0.2 for oe^URz and X a = 0.05 for a € 7£ 3 . 



0/ t/iem returned back to their own villages. This is, by the way, essentially what our 
results show. 



In particular, in Figures |2]j4] the parameter \ a is taken to be equal, A a = 0.05, in all 
the cells of 1Z. On the other hand, in Figures [5]{7j, \ a = 0.05 in 1Z 3 while \ a = 0.2 in 7Zi 
and 1Z 2 . This difference is useful to model the fact that Si and S 2 most probably interact 
where they live, rather than on the way. 

All these figures share a common feature: they all show that Si leaves IZi, moving 
towards 7Z 2 , while only a small part of S 2 moves towards TZ\. This is related to the value 
of the parameters 7 a and 75, as well as the p a ^ which were fixed at the very beginning, 



accounting for the diffusion in the model, see (3.3). Since 7 a > %, it is clear that Si has 



a larger mobility than S 2 . This is exactly what all the figures show. Figures 2(b) and 
3(b) show that, when the density of S\ in 7Z 2 approaches that of S 2 , S 2 reacts very fast 
in two ways: their birth rate increases very fast (since its density increases), and they 
start rejecting somehow the members of Si (since the density of Si decreases). After this 
first reaction, we see that, from time to time, a certain amount of people of Si goes back 
to TZi (presumably, after reaching some well-being). We see that in TZ 2 the density of S 2 
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Figure 7: Evolution of local densities. Africa (solid line for Si and dashed line for S2). Africa: top-left 
(a); Europe: top-right (b); Mediterranean Sea: bottom (c). j a — 0.1, 76 = 0.004, w£ = 1, = 2 
Va € ft; A Q = 0.2 for a e fti U ft 2 and A Q = 0.05 for a G ft 3 . 
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Figure 8: Evolution of local densities (solid line for Si and dashed line for 52). Inside 1Z C (left), 
Outside TZ C (right). j a = 0.1, 7 b = 0.004, X a = 0.2 for a £ R. First row: w£ = 1, w b a = 0.3, Va € ft. 
Second row: w£ = 1, w£ = 1, Va € ft. Third row: lo* = I, uj b a = 3, Va € ft. 
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stays almost always larger than that of Si, while in TZi the density of S 2 is always very 
low: rich people do not go in the poor area! Moreover, a lot of people of both populations 
are in TZ^: they travel, not necessarily moving from TZi to 1Z 2 or vice-versa. Incidentally 
we observe that, because of our interpretation in terms of richness of populations, the 
parameters 7 a and 7&, which are directly proportional to the mobility of the species, can 
also be seen as inversely proportional to their richness: the larger the value of j a , the 
poorer the species, and, consequently, the larger the will to go away from the related cell! 

The plots also suggest that the w's measure a sort of inertia of the two populations: 
increasing the value of, say, produces a less oscillatory behavior of S 2 , as we can see 
from Figures [2]-[4j Analogously, we have checked that increasing the value of u a , produces 
a more static behavior of Si. We should also mention that our numerical computations 
for L > 11 confirm our conclusions, showing that the size of 1Z is not important, except 
for the time needed to move from TZi to TZ 2 , as it is natural. Figures [5]-[7] show much faster 
oscillations in the densities of Si and S 2 than those in Figures [2]-[4| in particular in the 
regions TZi and 7Z 2 . This is due to the fact that, in these regions, the interaction parame- 
ters between the populations, A a , are taken larger than before. Hence, their densities can 
change faster than before: the interaction between Si and S 2 is more important than the 
diffusion of the populations! 

It should also be stressed that the somehow oscillatory behavior which is observed in 
many figures (and which can also be tested taking a larger time interval), reflects what 
observed in [20J in a different context, and is, in a sense, intrinsically related to the fact 
that we are dealing with a closed physical system ruled by linear ordinary differential 
equations possessing quasi-periodic solutions. The way our approach can be extended to 
include real damping effects is widely discussed in [2T], and is based on the introduction 
of a suitable reservoir interacting with the original system. In this case it is possible to 
check that the time dependence of the number operator a* a (t)a a (t) can be written as the 
product of an (essentially) oscillating self-adjoint operator xJ l (t)x Q ,(t), times a decaying 
factor: 

27T7 2 

a ] a {t)a a {t) = e'^ 1 x ] a {t)x a {t), 

where Q is a parameter of the hamiltonian of the reservoir, see [21] for a full analysis in 
a different context. 
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4 Competition between populations 



The same hamiltonian H introduced in Section 3 can be used in the description of com- 
petition between two populations Si and <S 2 . The equations for the two populations are 
again those in (3.4). The values of the p a ^ are chosen exactly as in the previous migration 
model: a component of S\ or 5 2 can move from one cell to a neightboring cell. Hence, 
the equation of motion can be written as in (3.10), and the solution is given in (3.11). 
The difference consists here in the choice of the parameters and of the initial conditions. 
In fact, in this case, we are no longer necessarily interested in having, at t — 0, the two 
populations localized in different regions of 1Z. Therefore, we consider here again a square 
lattice, with L = 11, in which both «Si and 5 2 are localized in a central region 7Z C of three 
by three cells, so that they are forced to interact between them from the very beginning. 

In Figure [8| in each row, we plot the local densities of Si (solid line) and <S 2 (dashed 
line) inside (left) and outside (right) 1Z C . Different rows correspond to different values of 
u b a , while all the other parameters coincide. We have chosen two significantly different 
values of 7 a and 7& to give the two populations different mobilities: since 7 a = 0.1 ^> 7& = 
0.004, Si is expected to move much faster than 5 2 , and this is exactly what we observe 
in the figure. Moreover, we have already seen that i>j a a and u h a play the role of inertia of 
the populations in the different cells. Hence, we expect that the higher is the ratio -f, 
the smaller will be the reaction of Si compared with that of iS 2 . These features are all 
evident in Figure [8] Si tends to move away from 1Z C faster (or even much faster) than 5 2 . 
Moreover, going from the first row {u^ = 1, uj b a = 0.3, Va G 1Z) to the last one (cu^ = 1, 
u b a — 3, Va £ 1Z), it is clear that the tendency of 5 2 to move away from 1Z C decreases 
more and more, even if its individuals keep on moving along 1Z. 

Particularly interesting is the second row where the density of S2 in 7Z C first increases 
very fast, while that of Si decreases: this can be considered as the evidence of a bigger 
efficiency of 5 2 compared with that of <Si, which is forced by 5 2 to leave 7Z C . For instance, 
thinking of Si as preys and of 5 2 as predators, we can say that the preys run very fast away 
from the region where the predators are localized. Hence 7 a and 76 can be considered, 
other than diffusion coefficients, as a sort of inverse ability of the two populations: since 
7 6 _1 ^> 7~ 1 , S2 is much stronger than Si, and the preys are killed significantly by the 
predators or, if they survive, run away from TZ C . 

Again, due to the absence of a reservoir, it is not expected any decay for large t, and 
this is exactly what the plots show. In order to have such a decay, the reservoir must be 
considered inside the model. This is reasonable, since such a reservoir can play the role 
of all the interactions that Si and 5 2 may experience other than the mutual interaction 
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(lack of food, other predators, cold winters, hot summers, . . . ). 

5 Conclusions 

In this paper, we have used fermionic operators to describe the dynamical behavior of 
two populations of individuals subjected to a certain diffusion. The model considered 
here is based on a quadratic hamiltonian, so that the resulting equations of motion are 
linear and the densities of the populations can be deduced analytically. We have used 
this rather general hamiltonian in two different contexts: first, we have described the 
dynamics of migration of two populations originally spatially separated. In this case we 
have seen that the migrants, which are originally well localized in a (poor) region of our 
lattice, move towards richer zones. This movement is driven not only from the general 
form of the hamiltonian but also by the choice of the parameters of the model, which 
therefore acquire a precise meaning. The same hamiltonian, with different choices of the 
parameters and of the initial conditions, has also been used in the description of the 
competition between two populations, like in a predator-prey simple system, and we have 
shown that again a reasonable and interesting dynamics can be deduced. 

Compared with other approaches and results, it is worth to be underlined that the use 
of fermionic operators automatically ensures the coexistence of the competing populations 
in the same environment. 

We are aware that the analysis here considered can be seen as a first step towards the 
construction of a more complete and satisfactory model of interaction between popula- 
tions. For instance, damping and/or nonlinear effects should be taken in consideration. 
Also, other possible topologies of the region 1Z may give interesting results. These are 
just part of our plans for the future. 

Appendix: Few results on the number representation 

We discuss here few important facts in quantum mechanics and in the so-called number 
representation, paying not much attention to mathematical problems arising from the 
fact that the operators involved might be unbounded, since this class of operators is 
not relevant for the applications proposed in this paper. More details can be found, for 
instance, in (27J EE], as well as in [ISl [IHl H3 HSl HH 1201 121] • 

Let Ti be an Hilbert space, and B(7-L) the set of all the bounded operators on 'K. 
Let S be our physical system, and A the set of all the operators useful for a complete 
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description of S, which includes the observables of S. For simplicity, it is convenient to 
assume that A coincides with B{H) itself. The description of the time evolution of S 
is related to a self-adjoint operator H = W which is called the Hamiltonian of S, and 
which in standard quantum mechanics represents the energy of S. We will adopt here the 
so-called Heisenberg representation, in which the time evolution of an observable X G A 
is given by 

X(f) =exp(iHt)Xex-p(-iHt), (A.l) 
or, equivalently, by the solution of the differential equation 
dX(t) 



dt 



iexp(iHt)[H,X]exp(-iHt) = i[H,X(t)\, (A.2) 



where [A, B] := AB — BA is the commutator between A and B. The time evolution 
defined in this way is a one-parameter group of automorphisms of A. 

An operator Z e A is a constant of motion if it commutes with H . Indeed, in this 



case, equation (A.2) implies that Z(t) = 0, so that Z(t) = Z for all t. 

In our previous papers [T5l [16], [TFJ [18], HHJ [20J, EJJ , a special role was played by the so- 
called canonical commutation relations. Here, these are replaced by the so-called canonical 
anti-commutation relations (CAR): we say that a set of operators {a^, a\,£ — 1, 2, . . . , L} 
satisfy the CAR if the conditions 

{at, al} = 6 in l, {ai,a n } = {a\,a ] n } = (A.3) 

hold true for all £, n — 1, 2, . . . , L. Here, 1 is the identity operator and {x, y} := xy + yx 
is the anticommutator of x and y. These operators, which are widely analyzed in any 
textbook about quantum mechanics (see, for instance, [271 EH]) are those which are used 
to describe L different modes of fermions. From these operators we can construct fii = a\a,£ 
and = Y2e=i ^ii which are both self-adjoint. In particular, fi£ is the number operator for 
the £-th mode, while N is the number operator of S. Compared with bosonic operators, 
the operators introduced here satisfy a very important feature: if we try to square them 



(or to rise to higher powers), we simply get zero: for instance, from (A.3), we have aj = 0. 
This is related to the fact that fermions satisfy the Fermi exclusion principle |28j . 

The Hilbert space of our system is constructed as follows: we introduce the vacuum of 
the theory, that is a vector ipo which is annihilated by all the operators af. a^o — for 
all I = 1, 2, . . . , L. Then we act on ipo with the operators a\ (but not with higher powers, 
since these powers are simply zero!): 

<fim,n 2 ,..,n L ■= (4r (4r • • • (4)"^o, (A.4) 
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n,£ = 0,1 for all £. These vectors form an orthonormal set and are eigenstates of both 
and N: h^m,n 2 ,...,n h = n^ ni ,n 2 ,...,n h and N(p nu n 2 ,...,n L = N(p ni>n2t ,„ )nL , where N = Y,t=i n t- 
Moreover, using the CAR, we deduce that hi (a^ ni)n2j ... jnz J = (ne — l)(aeip ni ,n 2 ,...,n L ) and 
fii (ci\(p ni)n2j __, >n[ /\ = ( n £ + 1) { a \ <fni ,n 2 ,...,n L ) j f° r all £. The interpretation does not differ 
from that for bosons, [T3], and then and a\ are again called the annihilation and the 
creation operators. However, in some sense, a\ is also an annihilation operator since, 
acting on a state with = 1, we destroy that state. 

The Hilbert space % is obtained by taking the linear span of all these vectors. Of 
course, % has a finite dimension. In particular, for just one mode of fermions, dimi^K) = 2. 
This also implies that, contrarily to what happens for bosons, the fermionic operators are 
bounded. 

The vector ip ni ,n 2 ,...,n L in 

as 

^ni,n 2 ,...,n L (X) = (if ), (A.5) 

where ( , ) is the scalar product in Ti. As we have discussed in [151 [13 HO HH1 UHl 1201 12T] . 
these states are used to project from quantum to classical dynamics and to fix the initial 
conditions of the considered system. 



A. 4) defines a vector (or number) state over the algebra A 
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